/*
 * Copyright (c) 1997, 1998, 2003, 2006 Aleksandar Samardzic
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are
 * met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. The name of the author may not be used to endorse or promote products
 *    derived from this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT,
 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */

#include "sphere.h"
#include "data_access_structure.h"
#include "voxel_grid.h"

/* Sphere_Intersect - find intersection of sphere with a line */
Bool
Sphere_Intersect(this, pLine)
     PSphere        *this;
     PLine          *pLine;
{
	double         *p0 = pLine->p0,
	    *dp = pLine->dp;
	double          radius = this->radius,
	    *center = this->center;
	PVector         diff;
	double          a,
	                b,
	                c;
	double          discr,
	                rec,
	                ext;
	double          t;

	Vector_Sub(diff, p0, center);
	/* relation |t*dp-diff|=r holds, which extends to square equation
	 * for t */

	/* calculate coefficients of equation */
	a = dp[X] * dp[X] + dp[Y] * dp[Y] + dp[Z] * dp[Z];
	b = 2 * Vector_DotProduct(diff, dp);
	c = diff[X] * diff[X] + diff[Y] * diff[Y] + diff[Z] * diff[Z] -
	    radius * radius;

	/* calculate root components */
	rec = 1 / (2 * a), ext = -b * rec;
	discr = ext * ext - 2 * c * rec;

	if (discr < 0)		/* imaginary roots */
		return FALSE;

	/* check if t is in range [EPSILON,pLine->t-EPSILON) */
	discr = sqrt(discr);
	if ((t = ext - discr) < EPSILON)
		t = ext + discr;
	if (t < EPSILON || t >= pLine->t - EPSILON)
		return FALSE;

	pLine->t = t;
	return TRUE;
}

/* Sphere_CalcNormal - calculate sphere normal at specified point */
void
Sphere_CalcNormal(this, point, normal)
     PSphere        *this;
     PVector         point;
     PVector         normal;
{
	double          mul,
	                den;

	Vector_Sub(normal, point, this->center);
	Vector_Normalize(normal, mul, den);
}

/* Sphere_CalcBoundingVolume - calculate bounding volume for sphere */
void
Sphere_CalcBoundingVolume(this, pNode, pSlabs)
     PSphere        *this;
     PBoundingVolume *pNode;
     PSlabs         *pSlabs;
{
	Byte            slab,
	                numSlabs = pSlabs->numSlabs;
	double          (*d)[2] = pNode->d;
	double          dist;

	for (slab = 0; slab < numSlabs; slab++) {
		dist =
		    Vector_DotProduct(this->center,
				      pSlabs->pVectors[slab]);
		d[slab][LO] = dist - this->radius;
		d[slab][HI] = dist + this->radius;
	}
}

/* Sphere_Voxelize - voxelization algorithm for sphere according to
 * [Andr97a] */
void
Sphere_Voxelize(this, pGrid)
     PSphere        *this;
     PVoxelGrid     *pGrid;
{
	PVector         center;
	double          radius = this->radius;
	Byte            slab;
	double          _hCell = 1 / pGrid->hCell;
	double          A,
	                B,
	                A1,
	                a1,
	                B1,
	                b1;
	double          z0,
	                z_max,
	                z_incr;
	double          R;
	double          t;

	Vector_Assign(center, this->center);
	for (slab = X; slab <= Z; slab++)
		center[slab] -= pGrid->d[slab][LO];
	Vector_MulAssign(center, _hCell);
	radius *= _hCell;

	A = (radius >=
	     sqrt(3) / 2) ? (radius - sqrt(3) / 2) * (radius -
						      sqrt(3) / 2) : 0;
	B = (radius + sqrt(3) / 2) * (radius + sqrt(3) / 2);
	z_incr = radius / (((int) radius) + 1);
	z0 = center[Z], z_max = center[Z] + radius;
	for (center[Z] -= radius; center[Z] <= z_max; center[Z] += z_incr) {
		A1 = A - ((center[Z] - z0) * (center[Z] - z0));
		B1 = B - ((center[Z] - z0) * (center[Z] - z0));
		if (A1 < 0)
			A1 = 0;
		if (B1 >= 0) {
			a1 = sqrt(A1), b1 = sqrt(B1);
			for (R = a1; R < b1; R += 1)
				VoxelGrid_EnumerateCellsCircleXY(pGrid,
								 center, R,
								 DBL_MAX,
								 &t,
								 VoxelGrid_Insert,
								 this);
			VoxelGrid_EnumerateCellsCircleXY(pGrid, center,
							 R - 1,
							 (B1 - A1) / 2, &t,
							 VoxelGrid_Insert,
							 this);
		}
	}
}
